On the polar caps of the Three Musketeers 



A. De Luca, P. A. Caraveo, S. Mereghetti, M. Negroni^, G.F. Bignami 

Istituto di Astrofisica Spaziale e Fisica Cosmica, 
Sezione di Milano "G.Occhialini" - INAF 
; v.Bassini 15, 1-20133 Milano, Italy 

O \ 

O ! delucaSmi . iasf . cnr . it 

(N 

o 

CD 

Q 



1,2 



0^ 



ABSTRACT 

XMM-Newton observations of PSR B0656+14, PSR B1055-52 and Geminga 
have substantially increased the statistics available for these three isolated neu- 
^ ' tron stars, so apparently similar to deserve the nickname of "Three Musketeers" 

^ ■ (Becker & Triimper, 1997). Here we shall take advantage of the EPIC statistics 

CSJ ' to perform phase resolved spectroscopy for all three objects. The phase-averaged 



spectrum of the three musketeers is best described by a three component model. 
O ■ This includes two blackbody components, a cooler one, possibly originating from 

• the bulk of the star surface, and a hotter one, coming from a smaller portion 

of the star surface (a "hot spot"), plus a power law. The relative contributions 



2 ■ of the three components are seen to vary as a function of phase, as the stars' 

^ ■ rotation bring into view different emitting regions. The hot spots, which have 

very different apparent dimensions (in spite of the similarity of the three neutron 
stars polar cap radii) are responsible for the bulk of the phase variation. The am- 
plitude of the observed phase modulation is also markedly different for the three 
sources. Another striking aspect of our phase-resolved phenomenology is the ap- 
parent lack of any common phase alignement between the observed modulation 
patterns for the two blackbody components. They are seen to vary in phase in 
the case of PSR B1055-52, but in anti-phase in the case of PSR B0656-I-14. These 
findings do not support standard and simplistic models of neutron star magnetic 
field configuration and surface temperature distribution. 
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PSR B1055-52) — stars: neutron — x-ray:stars 



1. Introduction 

PSR B0656+14, PSR B1055-52 and Geminga are isolated neutron stars showing similar 
periods and period derivative. Hence their characteristic ages (few 100 thousand years), 
their inferred magnetic fields (few 10^^ G) and their energetics (-E^ot ~ 10^^ erg s~^) are 
very similar. Moreover, they are all nearby, making it easier to detect them at different 
wavelengths. Parallax distance measurements are available for Geminga (157 pc, Caravco 
et al., 1996) and PSR B0656+14 (288 pc, Brisken et al., 2003), while for PSR B1055-52 
the current distance estimate, based on the dispersion measure, is ~750 pc (Kramer et al., 
2003). In X-rays they shine both owing to thermal and non-thermal mechanisms. Black- 
body emission, from the majority of their surface, is seen to be complemented at higher 
energies by power law "tails" . They have faint optical counterparts, with visual magnitude 
m^ ~25.5 for Geminga (Bignami et al, 1987), ~25 for PSR B0656+14 (Caraveo et al., 1994) 
and ~24.9 for PSR B1055-52 (Mignani et al., 1997). Multicolor photometry, available for 
Geminga and PSR B0656-I-14, shows that their NUV magnitudes are broadly compatible 
with the extrapolations of their X-ray blackbody emissions, while a power law-like excess 
is apparent in the optical/NIR (Mignani et al., 2004 and references therein). While PSR 
B1055-52 and Geminga arc bona fide gamma-ray pulsars (Fierro et al., 1993; Bertsch et al., 
1992), for PSR B0656H-14 a tentative detection awaits confirmation (Ramanamurthy et al., 
1996). 

Here we shall concentrate on the analysis of the XMM-Newton data collected both in 
imaging and timing mode. Such exceptional harvest of X-ray photons over the ample energy 
interval covered by the EPIC cameras allows us to use a new tool for the study of the X-ray 
behaviour of the "Three musketeers" : phase resolved spectroscopy. 

While we have already published the phase-resolved analysis of the EPIC/pn data on 
Geminga (Caravco et al. 2004), here we shall include also this source, in order to underline 
similarities and differences in the phase resolved X-ray behaviour of three otherwise very 
similar objects. 



^ Based on observations with XMM-Newton, an ESA science mission with instruments and contributions 
directly funded by ESA member states and the USA (NASA). 
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1.1. Review of the X-ray observations of the Three Musketeers 

The EPIC data presented below come after two decades of X-ray observations of isolated 
neutron stars, starting with their discovery as soft, probably thermal sources by the Einstein 
Observatory. 

1.1.1. Geminga 

IE 0630+178 was discovered by the Einstein satellite while covering the gamma-ray error 
box of the unidentified source Geminga (Bignami et al., 1983). The discovery of a common 
237 msec pulsation clicked the identification, the first ever, of an unidentified gamma-ray 
source with its X-ray counterpart (Halpern & Holt, 1992; Bertsch et al., 1992). The first 
good quality spectral data on Geminga were collected by ROSAT. Halpern and Ruderman 
(1993), using 7,911 photons in the 0.07-1.50 keV energy band, concluded that the spectrum 
was well described by two black-body curves at temperatures of 0.5 x 10^ K and 3x 10^ K. 
While the cooler component appeared to come from the full surface of the neutron star, the 
hotter one could have been coming from a spot covering just 3x10^^ of the neutron star 
surface, possibly due to a heated polar cap. A second, longer ROSAT observation, yielding 
18,500 soft photons, was analysed together with an observation by ASCA, providing the 
first 1,750 "hard" (> 2 keV) counts. Such a combination seemed to favour a different, 
composite interpretation, where, a power-law, non-thermal emission was added to the low 
temperature black-body radiation (Halpern and Wang, 1997). The power-law photon index 
was still poorly constrained, ranging from 2.2 to 1.5. Using a longer ASGA observation 
with 4,800 photons, Jackson et al. (2002) confirmed the composite nature of Geminga's 
spectrum and refined the power-law photon index value to 1.72±0.1. They tried also to 
perform some phase-resolved spectroscopy, but could only use two phase intervals. They 
found the spectrum of the interval including the peak to be slightly softer than the rest of 
the light curve, but concluded that "the difference is only marginally significant" . 

1.1.2. PSR B0656+U 

The definitely brighter PSR 0656+14 was first detected by the Einstein satellite (Cor- 
dova et al., 1989). A ROSAT PSPC observation allowed Finley et al. (1992) to detect the 
X-ray pulsation and to measure a pulsed fraction of 14% ± 2%. Further pontings were 
then carried out with the ROSAT PSPC detector in 1992 (Oegelman, 1995), for a total 
exposure time of about 17 ksec, collecting ~32,000 photons in the 0.1-2.4 keV band. The 
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overall ROSAT dataset was analyzed by Possenti et al. (1996), who found the bulk of the 
emission to be of thermal origin, well described by a blackbody curve (T~ 9 x 10^ K) orig- 
inating from a large part the star surface (emitting radius ~14 km, assumig a distance of 
760 pc, corresponding to ~5.3 km at the parallactic distance of 288 pc). A second spectral 
component was required to describe the higher energy part of the spectrum, as well as to 
explain a significant change of the pulse profile with energy It was not possible however 
to discriminate between a second blackbody component (T~ 1.9 x 10^ K) originating from 
a region a few hundred times smaller and a steep (F ~ 4.5) power law. Greiveldinger et 
al.(1996) coupled the ROSAT dataset to an ASCA observation which yielded ~2,000 pho- 
tons in the 0.5-10 keV band. Their best fit used a three component model: the sum of two 
blackbodies (best fit parameters very similar to the results of Possenti et al. 1996) and of 
a power law with photon index r=1.5±l.l. More recently, PSR B0656+14 was observed 
with Chandra, both with ACIS and with the LETG. The ACIS data (~45,000 photons in 
the 0.2-6 keV range) confirmed a three-component model to yield the best description of the 
spectrum (Pavlov et al., 2002), consistent of the sum of two blackbodies (Ti ~0. 85x10^ K, 
Ri ~12 km; T2 ~1.65xl0^ K, R2 ~1 km, using a distance of 500 pc; assuming the distance 
of 288 pc, Ri ~7 km and R2 ~0.6 km). The high spectral resolution of the Gratings allowed 
Marshall & Schultz (2002) to exclude the presence of significant features superimposed on 
the thermal continuum in the softer band (0.15-1 kcV). The parameters of their best fitting 
two-blackbody model are Ti -^0.8x10^ K, Ri ~22.5 km; T2 ~1.5xl0^ K, R2 ~1.7 km, 
assuming a distance of 760 pc, corresponding to Ri ~8.5 km and R2 ~650 m at the distance 
of 288 pc. 



1.1.3. PSR B1055-52 

After the Einstein observatory discovery of X-ray emission from this radio pulsar (Cheng 
& Helfand, 1983), PSR B1055-52 was observed with ROSAT, both with the HRl (8.6 ksec 
yielding ~ 570 source photons) and with the PSPC (15.6 ksec, for a total of ~ 5500 source 
photons) in 1990-1992 (Ocgelman & Finley, 1993). The timing analysis unveiled the source 
pulsation, with a pulsed fraction increasing from ~11% for energies <0.5 keV to ~63% above 
0.5 keV. The spectrum was best described by a two component model. A blackbody with 
temperature of ~ 8 x 10^ K, coming from a large portion of the surface, accounts for the bulk 
of the X-ray emission; a second component was required, either a second, hotter (T~ 3.6 x 10^ 
K) blackbody coming from an area a few 10^ times smaller, or a steep (F ~ 4) power law. 
An ASCA observation could add only ~200 photons in the 0.5-10 keV range (Greiveldinger 
et al., 1996). A Chandra ACIS observation was taken in 2000 (42 ksec exposure); results 
were summarized by Pavlov et al. (2002). A simultaneous fit to Chandra and ROSAT data 
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required a three component model, consisting of the sum of two blackbodies (Ti ~8.3xl0^ 
K, Ri ~12 km; kT 2 ~ 1-6 x 10^ K, R2 ~800 m assuming the distance to be 1 kpc; Ri ~9 km, 
R2 ~600 m at the revised distance of 750 pc) and of a power law (F ~1.7). The above authors 
reported the presence of a significant variation of the hot blackbody component, observed 
through phase-resolved spectroscopy, but they did not present detailed results. A similar 
three component model has been used by Becker & Aschenbach (2002) in their analysis of 
XMM-Newton data. Spectroscopy was performed using MOS data only, which were fitted 
together with ROSAT data. Their best fit was obtained with Ti ~7.1xl0^ K, Ri ~31 km, 
T 2 ~ 1.4 X 10^ K, R2 ~2.6 km assuming the distance to be 1 kpc (at the revised distance of 
750 pc the emitting radii are Ri ~15.5 km, R2 ~1.3 km); the power law photon index was 
found to be 1.9±0.2 and the interstellar absorption Nh ~ 2.3 x 10^° cm~^. 

2. The XMM-Newton/EPIC data sets 

Geminga, PSR B0656+14 and PSR B1055-52 were observed by XMM-Newton as Guar- 
anteed Time targets. The complete journal of observations is reported in Table 1. The data 
are now available in the XMM-Newton Science Archive. 

While the MOS cameras (Turner et al., 2001) were operated in Pull Frame mode in order 
to image on the full field of view of the telescopes (15' radius), the pn detector (Striider et al., 
2001) was used to time-tag the photons, either in Small Window mode (6 ms time resolution, 
imaging on a 4' x 4' field) or in Fast Timing Mode (0.03 ms resolution, no spatial information 
along detector columns). Unfortunately, as pointed out by Becker & Aschenbach (2002), the 
use of the potentially promising Fast Timing mode produces, as a byproduct, a background 
significantly higher than in Small Window mode, owing to the peculiar readout mode and to 
the collapse of data along the CCD columns. Such an abnormally high background lowers 
the signal-to-noise above a few keV. 

2.1. Data reduction and background screening 

All the data reduction was performed using the most recent release of XMM-Newton 
Science Analysis Software (SASv6.0.0). The raw Observation Data Files (ODFs) were pro- 
cessed using standard pipeline tasks {epproc for pn, emproc for MOS data). We selected 
events with PATTERN 0-4 and PATTERN 0-12 for the pn and the MOS, respectively. 

Particular care was devoted to reduce the instrumental background in the linearized 
event lists. First, an accurate screening for soft proton fiare events was done, following 
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the prescription of De Luca & Molendi (2004). We computed Good Time Intervals (GTIs) 
setting a threshold of So" from the quiescent rate on the 0.5-12 keV count rate for both 
the pn and the MOS detectors. A more stringent threshold was adopted in the case of pn 
Fast Timing mode observations, since the collapse of data along the CCD readout direction 
increases (by a factor ~ 30) the background count rate in the source extraction region. 

The pn Fast Timing mode data required a further cleaning. Such operating mode is 
affected by a peculiar flaring background component ("Soft Flares", see Burwitz ct al., 2004) 
ultimately due to the interaction of charged particles (possibly heavy ions) with the CCD. 
Such interactions produce short (0.1-0.5 s time scale) and very intense bursts of events, 
with typical energies of ~ 0.22 keV for singles (mono-pixel events) and of ~ 0.45 keV 
for doubles (bi-pixel events). To remove such background noise, which would hamper the 
study of the low-energy emission of our targets, ad-hoc GTI files for singles and doubles 
were created, extracting 1-sec binned light curves in the 0.2-0.3 keV and 0.4-0.5 keV energy 
ranges, respectively, and then using the same algorithm adopted for the soft proton case. 

The overall GTI filter for each pn Fast Timing observation was then obtained by in- 
tersecting the soft proton GTI, the soft flares GTI for singles and the soft flares GTI for 
doubles. The net exposure times (after dead-time correction) for the cleaned event hsts are 
reported in Table 1. In the case of PSR B1055-52 the observation was splitted between two 
XMM revolutions (186 and 187). Event lists obtained with the same instrument in the same 
mode in different orbits were then merged using the SAS task merge. 

2.2. Source and background events selection 

To extract the source photons from the dataset taken in imaging modes (Full Frame 
for MOS and Small Window for pn) we selected a circle of 45" radius, containing ~ 90% of 
the counts; in Fast Timing mode we used a 17 pixel wide strip (4.1" pixel size), containing 
~ 85% of the source counts. Background photons were extracted from suitable regions on 
the same CCD chip containing the source: for Full Frame MOS images we need an annulus of 
2' inner radius and 4' outer radius; for the pn detector operated in small window we selected 
two rectangular regions oriented along the CCD readout direction and covering ~2 square 
arcmin; for pn used in Fast Timing mode we selected two stripes (17 and 12 pixels wide) 
away from the source region. 

The 0.25-8.0 keV count rate observed in MOSl camera for PSR B0656-M4 is of 0.9 
counts s~^. A modest pile-up is expected, owing to the slow CCD readout in the Full Frame 
mode (2.6 sec frame time). An analysis of the event PATTERN distribution, performed with 
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the SAS task epatplot, clearly showed such effect as an excess of bi-pixel events above 0.7 
keV. To solve the problem, we simply excluded the PSF core (inner 5", containing about 
35% of the source photons), where the pile up could be significant. No hints for pile-up in 
the resulting photon list was apparent in the corresponding event PATTERN distribution. 

In the case of PSR B1055-52, EPIC MOS images show a faint source at an angular 
distance of 32" in the NE direction (see Figure 17 of Becker & Aschenbach, 2002). Such 
source has a hard spectrum (possibly a background AGN); its flux, negligible wrt. the 
pulsar below 2 keV, becomes comparable to the pulsar's one above 3 keV. In order to avoid 
possible contaminations in the study of the pulsar high energy emission wc decided to exclude 
such source from the pulsar extraction region, using a 10" arcsec radius circle (containing 
~ 60% of the counts at 5 keV) as a geometrical mask. Such a selection was not possible in 
pn Timing mode, since along the RAWX direction (perpendicular to the readout direction), 
where spatial information is mantained, the angular separation of the two sources is <4 
pixels. Therefore, the flux measured by the pn above ~3 keV is possibly contaminated by 
some contribution from such background source and its absolute value should be taken with 
caution. 

Background-subtracted count rates for the three neutron stars are reported in Table 1, 
together with the number of collected photons. 

3. EPIC Data analysis 

With a number of time-tagged photons (see Table 1) more than doubling all previous 
statistics, EPIC offers now the first chance of meaningful phase-resolved spectroscopy for the 
three musketeers. 

Following the procedure outhned by Caraveo et al. (2004): 

• first, we address the phase averaged spectral shape of the three objects to obtain useful 
pieces of information on the components needed to fit their overall spectra; 

• next, we perform the timing analysis to search for the best period and to perform the 
phase alignement; 

• then, we divide the folded light-curves in 10 phase intervals and extract the spectra 
corresponding to each phase interval; 

• finally, each spectrum is fitted using the same components found to best describe 
the integrated spectrum and the corresponding phase-resolved spectral parameters are 
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computed. 

In the following section we will discuss in detail each step for the cases of PSR B0656+14 
and PSR B1055-52. The case of Geminga was described by Caraveo et al.(2004); we will 
report here the results to ease a synoptic view of the phenomenology of the three musketeers. 

3.1. Phase-integrated spectral analysis 

Spectra for the source and background events, extracted using the regions described 
in Sect. 2.2, were rebinned in order to have at least 40 counts/channel. We added a 5% 
systematic error to each spectral bin, to account for calibration uncertainties among different 
instruments/modes. Ad hoc redistribution matrices and effective area files were generated 
using the rmfgen and arfgen tasks of the SAS. The spectral analysis was performed with 
XSPEC vll.2. We selected 0.25 keV as lower bound for the spectral study, since below such 
an energy the calibration is still uncertain; the upper bound was selected on the basis of the 
observed signal-to-noise, which depends both on the sources' spectra and on the operating 
mode used. As discussed in Sect. 2.1, the spectra collected by the pn in Fast Timing mode 
suffer from a rather high background, which hampers the study of our targets at high energy, 
where such sources are rather faint. For PSR B0656+14 spectra obtained in Fast Timing 
mode turned out to be useful up to 2 kcV, while for PSR B1055-52 we could include data up 
to 6 keV. Spectra from MOSl/2 Full Frame mode, as well as from pn sw mode, were used 
up to 8 keV. MOS and pn spectra were fitted simultaneously, leaving a cross-normalization 
factor as the only free parameter. 

No significant spectral features are detected superimposed on the continuum of the three 
musketeers, neither in emission nor in absorption. As in the case of Geminga (Garaveo et al., 
2004) (Fig.l), the best fitting model is found combining two blackbody curves and a power 
law. This confirms the findings of Pavlov et al. (2002), based on Ghandra observations of 
PSR B0656+14 and PSR B1055-52, and Becker & Aschenbach (2002), who analyzed XMM- 
Newton (MOS only) data for PSR B1055-52. Using the best fitting parameters reported in 
Table 3 we obtain x^=l-l (368 d.o.f.) for PSR B0656+14 (Fig.2) and x^=l-0 (327 d.o.f.) for 
PSR B1055-52 (Fig. 3). The use of a second power law instead of the hotter blackbody yields 
photon indices of ~ 6.5 for PSR B0656+14 and ~ 5.6 for PSR B1055-52, while lowering 
the fit quality {xl >l-5 and xl >1-1 for the two cases, respectively). We note that a 
single-temperature, magnetized neutron star atmosphere model (Zavlin et al. 1996) cannot 
adequately reproduce the profile of the low-energy part of the spectrum of PSR B0656-I-14 
and PSR B1055-52, similarly to the case of a single blackbody curve. Moreover, such model 
requires a very large emitting surface (with emitting radii of ~100 km - this is true also for 
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Geminga) . Therefore such model cannot provide a satisfactory description of the spectra of 
the three musketeers. 

A synoptic plot of the spectra of the three musketeers is shown in Fig. 4. It is easy to 
note the definitely higher signal to noise in the high energy portion (>2 keV) of the spectrum 
of Geminga, in spite of a flux only slightly higher in such energy range. This is mainly due 
to the different operating modes used for the pn detector in the observations of the three 
targets, especially to the much lower background in the pn small window mode wrt. Fast 
timing mode. The power law component in the spectra of PSR B0656+14 and PSR B1055- 
52 is constrained by the MOS spectra and (in the case of PSR B0656+14) by the short 
pn small window observation (see the caption to Fig. 4 for further details). We summarize 
the results of the EPIC spectral fits for the three musketeers in Table 3, where errors are 
computed at 90% confidence level for a single interesting parameter. The lower temperature 
blackbody (hereafter "cool blackbody") is associated to an emitting region with an area 
compatible with the full surface of the neutron star. The higher temperature component 
( "hot blackbody" ) is seen to originate from a much smaller area. To visualize the correlation 
between different spectral parameters we computed confidence contours (at 68%, 90% and 
99% levels for two parameters of interest) for the interstellar column density Njj vs. the cool 
blackbody emitting area; this is shown in the insets in Fig. 1, Fig. 2 and Fig. 3. 

As shown in Fig. 4, the overall shape of the spectrum of the three musketeers is very 
similar. However remarkable differences in the luminosity of different spectral components 
are immediately evident (see Table 3). In particular, we note that the luminosity of the hot 
blackbody in the case of Geminga is ~ 2 orders of magnitude lower than in the other two 
cases. At variance with PSR B0656+14 and PSR B1055-52, for Geminga the contribution of 
the power law dominates over the hot blackbody in the energy range where the hot blackbody 
has its maximum. 

The parameters best fitting the EPIC spectrum of PSR B1055-52 (see Table 3) are 
fully consistent with the results reported by Becker & Aschenbach (2002) and by Pavlov et 
al.(2002), who included ROSAT data in their spectral fits to EPIC/MOS and Chandra/ACIS 
data, respectively. The same is true in the case of Geminga: although we use a 3 component 
model, the results for the Cool blackbody component (temperature and emitting radius) are 
found to be very similar to the ROSAT ones (see e.g. Halpern & Ruderman, 1993). The 
case of PSR B0656+14 is slightly different: the interstellar column density (4.3 ± 0.2 x 10^° 
cm~^) resulting from the symultaneous fit to the pn/Fast Timing, pn/Small Window and 
MOSl/FuU Frame data is larger than the values (~ 2 x 10^° cm^^) obtained by Pavlov 
et al.(2002) and Marshall & Schulz (2002) (based on Chandra ACIS/LETG data) and by 
Possenti et al.(1996) (based on ROSAT data). Thus, our analysis of the EPIC data results 
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in a brighter underlying continuum at low energies, leading to a lower temperature coupled 
to a larger emitting surface for PSR B0656+14. Addressing such a discrepancy, possibly 
due to subtle cross-calibration issues, is beyond the scope of our fitting exercise, which was 
performed to provide the starting point for the phase- resolved spectroscopy of the EPIC 
data. 



3.2. Timing analysis 

As a first step, wc have studied the high time resolution pn data in order to derive the 
period of our targets. Source photons were selected from the regions described above, in 
the energy range 0.15-8 keV for pn Small Window mode and 0.15-2 keV for pn Fast Timing 
modc^. Photon time of arrivals were converted to the Solar System barycenter using the 
SAS task barycen and folded using 10 bins in a range of trial periods around the expected 
values. A very significant detection of the pulsation was obtained in each case. To determine 
with higher accuracy the period value and to evaluate the corresponding error, we followed 
the prescription of Leahy (1987). 

The results are reported in Table 2, where they are compared with the extrapolations of 
the radio ephemeris, selected to be the nearest to the XMM-Newton observation. The period 
for the radio-quiet Geminga is compared to the 7-ray ephemeris of Jackson et al. (2002), 
based on the analysis of the complete EGRET dataset. A perfect consistency between the 
expected and measured values is apparent. 

In order to align in phase the X-ray light curves with the radio (or gamma) ones, we 
have folded XMM-Newton data using the radio timing solutions for PSR B0656-I-14 and PSR 
B1055-52 and the gamma one for Geminga. The absolute accuracy of the XMM-Newton clock 

is ~ 500 /is, as stated by the calibration team (Kirsch, 2004) and by recent investigations 
(e.g. Becker ct al., 2004). The extrapolation of the radio ephemeris to the time of the XMM 
observations yields a phase uncertainty of ~0.01 for PSR B0656+14 (corresponding to 1/5 of 
phase bin in Fig. 5) and of ~0.003 (1/20 of phase bin in Fig. 5) for PSR B1055-52, for which 
the clock accuracy is the dominant factor, since the radio data are almost simultaneous to 
the X-ray ones. 

A larger uncertainty, ~ 0.15 in phase, is present in the case of Geminga, as discussed 
by Caraveo et al. (2004), owing to the much longer time separation between the EGRET 



^Time-tagging of photons is reliable over the whole detector sensitivity range; therefore we include in the 
timing analysis also low energy (E<0.25 keV) photons, not used for spectroscopy (see Sect. 3.1) 
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and XMM-Newton observations. 

The folded pulse profiles (background-subtracted) of our targets using the overall EPIC 
energy range are shown in Fig. 5. To ease their comparison, all the lightcurves have been 
plotted setting phase to the X-ray pulse maximum as seen in the overall energy band. The 
phases of the radio peaks (PSR B0656-I-14 and PSR B1055-52) and 7-ray peaks (Geminga and 
PSR B1055-52) are also shown. Pulsed fractions have been computed as the ratio between 
the counts above the minimum and the total number of counts, deriving the corresponding 
uncertainties from the propagation of the statistical errors in the folded light curves. 

We have then studied the variations of the pulse profiles with energy. The light curves 
relative to four energy intervals are shown in Fig. 6 (PSR B0656+14), Fig 7 (PSR B1055-52) 
and Fig. 8 (Geminga). Note that the energy ranges used are different for different targets, 
since for each source the energy intervals have been chosen in order to visualize the light curve 
relative to the cool blackbody, transition region, hot blackbody and power law. The lower 
panels of figures 6 and 7 show the radio profile of PSR B0656+14 and PSR B1055-52 at 0.67 
GHz^. For the case of Geminga (Fig. 8), the 7-ray light curve^ is shown. Ten phase intervals, 
selected a priori, to be used for phase-resolved spectroscopy (see Sect. 3.3) are marked in the 
upper panel. Our three neutron stars exhibit vastly different phenomenologies. The main 
results may be summarized as follows (see also the captions to Fig. 6, Fig 7 and Fig. 8): 

• For PSR B0656-I-14 a single pulse per period may be seen in each energy band. The 
pulse shape below 1.5 keV is broadly sinusoidal; a significant change in pulse profile 
occurs between the energy range dominated by the cool blackbody and the range 
dominated by the hot blackbody. As a consequence, in the intermediate range the 
pulsed fraction reaches its minumum value. At higher energy, the pulse profile seems 
sharper, although the low statistics does not allow to study its shape (nor the pulsed 
fraction) with great accuracy. Note that the light curve above 1.5 keV was obtained 
by folding data from the short pn Small Window mode observation; the other cases 
represent data from the longer Fast Timing observation. Gonsidering the overall hght 
curve (Fig. 5), the single radio peak lags the X-ray maximum by ~ 0.25 ±0.05 in phase 
(the uncertainty is dominated by the X-ray light curve bin width for the determination 
of the position of the X-ray peak), occurring very close to the X-ray minimum. In the 



•^Retrieved from the European Pulsar Network database, http://www.nipifr.bonn.mpg.de/div/pulsar/data/browser.html 

^this was obtained by selecting 9 EGRET Viewing periods in which the target was imaged close to the 
center of the field of view and having a good photon statistics; events with energy E>100 MeV were extracted 
from the source position; time of arrivals were converted to the solar system barycenter and folded with the 
ephemeris published by Jackson et al. (2002) 
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lower panel of Fig. 6 we mark the phases of the double-peaked pulse profile observed 
in the V band by Kern ct al.(2003). The position of the single 7-ray peak claimed by 
Ramanamurthy et al.(1996) is also indicated. 

• Also PSR B1055-52 shows a single pulse per period at all energies, with a markedly 
different shape in the softest band (below 0.35 keV) wrt. intermediate energy ranges 
(0.35-1.5 keV), while above 1.5 keV the low signal-to- noise hampers the study of the 
pulse profile. The pulsed fraction, which is seen to grow with energy, is remarkably high 
(~ 70%) in the band where the hot blackbody dominates the neutron star emission. 
The two peaks of the radio profile are seen to contour the X-ray peak at all energies; in 
the overall light curve the highest radio peak lags the X-ray one by ~ 0.2 ± 0.05 (the 
uncertainty is dominated by the X-ray light curve binning). The phase interval where 
the 7-ray pulse occurs lies between the X-ray maximum and the highest radio peak. 

• Geminga shows a different behaviour of the pulse shape as a function of energy. The 
single, broad peak seen at low energy, where the cool blackbody dominates, changes to 
a double peak towards higher energies, where the bulk of the emission is non-thermal 
(power law component). The pulsed fraction reaches a maximum of >50% in the 
energy range where the hot blackbody is more important. The significant uncertainty 
(~ 0.15) affecting the extrapolation of the EGRET ephemeris (see above) prevents any 
firm conclusion on the alignement between the X-ray and the 7-ray light curves. 

3.3. Phase-resolved spectroscopy 

The results of the spectral and timing analysis were then used to study the variation 
of the emission spectrum with pulse phase. For both PSR B0656-I-14 and PSR B1055-52 
we selected photons whose times of arrival fall into the 10 different phase intervals marked 
on Fig. 6 and Fig. 7. The time resolution of the pn Fast Timing mode (0.03 ms) and Small 
Window mode (5.6 ms) are adequate to perform such photon phase selection. We then 
extracted 10 phase-resolved spectra. For the case of PSR B0656+14 this was independently 
done for both Fast Timing data and Small Window data. Each spectrum was then rebinned 
in order to have at least 40 counts per bin. Since now we arc interested to study relative 
variations as a function of the pulsar phase, we didn't add systematic errors to the pn 
phase-resolved spectra. 

As already suggested by the analysis of the energy-resolved pulse profiles, the spec- 
tra of both sources do vary as a function of the pulsars' rotational phase. This is easily 
seen in the first column of Fig. 9a and 9b (PSR B0656+14), Fig. 10a and 10b (PSR 
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B1055-52), where the observed phase-resolved spectra are compared to the best fit model 
describing the phase-averaged ones and the deviations arc computed in units of statisti- 
cal standard errors (an animated version of Figures 9a/b, lOa/b and lla/b is available at 
http://www.mi.iasf.cnr.it/~deluca/3musk/). The variation is dramatic in the case of PSR 
B1055-52, owing to the high phase modulation of its emission (up to 70% in the 0.7-1.5 kcV 
range, see Fig. 7) while it is less spectacular (but still seen with high significance) for PSR 
B0656-I-14, due to its lower modulation (pulsed fraction <20% below 2 keV, see Fig. 6). Note 
that in the case of PSR B0656-I-14 phase-resolved spectra have been rebinned in order to 
ease the visibility of the modulation on the graphs' scale. PSR B0656-I-14 and PSR B1055-52 
seem therefore to share the same characteristics seen in Geminga, but their behaviours follow 
different trends (the case of Geminga is reported in Fig. 11a and lib). 

Turning now to the spectral fit of the phase-resolved spectra, no meaningful results could 
be obtained by allowing both the temperature and the normalization of the blackbodies to 
vary: a definite trend showing the maximum temperature in correspondence of minimum 
emitting area (and vice-versa) was systematically obtained in all cases. Such a modulation 
pattern is obviously driven by the strong correlation between the two spectral parameters, 
and cannot be trusted as reliable. 

Thus, wc have used the phase-integrated best fit model as a template to describe the 
phase-resolved data. The N// value, the temperatures of both blackbodies and the power 
law photon index were fixed to the values best fitting the phase-integrated spectrum. The 
phase modulation was then reproduced by allowing the normalization parameters of each 
spectral component to vary independently. For the case of PSR B0656-I-14, spectra from 
Small Window data and from Fast Timing data corresponding to the same phase interval 
were fitted simultaneously. 

As in the case of Geminga (Caraveo et al., 2004), such a simple approach (using only 
three free parameters) yielded a satisfactory description of the phase-resolved spectra also 
for PSR B0656+14 and PSR B1055-52. Reduced of ~ 0.8 - 1.4 were obtained in different 
phase intervals. 

The results of phase-resolved spectral fits are shown in right column of Fig. 9a and 9b, 
Fig. 10a and 10b (the case of Geminga is represented in Fig. 11a and lib): for each phase 
interval, the three best fitting spectral components are plotted, superimposed to the unfolded 
data points. Knowing the objects' distance values, we can then compute the emitting surface 
corresponding to each blackbody curve. To visuahze the evolution of the spectral parameters 
as a function of the pulsars' phase, we have plotted the blackbody radii, as well as the power 
law intensities, as a function of the objects' rotational phase. This is shown in Fig. 12, 
Fig. 13 and Fig. 14. 
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The behaviour of the three neutron stars as a function of their rotational phase is vastly 
different. The most important results may be summarized as follows: 

• The cool blackbody component shows a similar phase evolution for the three targets. 
The values of the blackbody radii follow a roughly sinusoidal profile with a modulation 
<10% (wrt. the average value) for PSR B0656+14 and for PSR B1055-52, of -15% 
for Geminga. 

• The hot blackbody components have a similar phase evolution with a sinusoidal profile. 
Hovever, we see striking differences both in the amplitude of their phase modulation 
and on their overall luminosity. While for PSR B0656+14 the modulation in the 
emitting radius wrt. the average value is <10%, similar to the value found for the cool 
blackbody component, in the case of PSR B1055-52 we see a 100% modulation, since 
the hot blackbody component is not seen in 4 out of 10 phase intervals. A similar, 
100% modulation is observed also in the case of Geminga, although in this case the hot 
blackbody component is seen to disappear in just one phase interval, and the profile 
of its phase evolution is markedly broader. As far as the radii values are concerned, 
we note that they span a factor of 30, ranging from the 60 m of Geminga to the ~2 
km of PSR B0656+14. 

• The study of the power law component is partially hampered, in the cases of PSR 
B0656+14 and PSR B1055-52, by the low signal-to-noise at energies above 1.5 keV, 
owing to the high background affecting the Fast Timing observations. It is therefore 
difficult to assess the shape of the pulse profile (single-peaked for PSR B0656+14?; 
double-peaked for PSR B1055-52?) and the actual modulation in such cases. Con- 
versely, for Geminga, thanks to the different pn mode, the power law component show 
a clear double-peaked profile, with a significant unpulsed fiux. 

• The relative phase evolution of the three spectral components is vastly different for the 
three neutron stars. This is particularly evident when looking at the two blackbody 
components. In the case of PSR B0656H-14 some sort of anti-correlation is observed, 
with the hot blackbody peaking in correspondence of the cool blackbody minimum. 
This is at odds with the behaviour observed for PSR B1055-52, where the two thermal 
components appear definitely correlated. Geminga shows an intermediate phenomenol- 
ogy: the cool and the hot blackbodies have a ~ 0.25 difference in their phase pattern. 

Of course, one could fit the phase-resolved spectra fixing the emitting radii to the average 
values and leaving the temperatures as free parameters. The goodness of the fits obtained 
following such an approach is comparable to that obtained previously fixing the temperatures 
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and the photon index. Not surprisingly, the evolution of the parameters, as a function of 
the phase, is very similar to the plots shown in Fig. 12, Fig. 13 and Fig. 14. However, such 
an approach renders less immediate the interpretation of the hot blackbody phase evolution 
especially in the case of PSR B1055-52, where the hot spot is not present for about half of 
the period. Maximum emitting radius translates into maximum temperature; absence of hot 
spot translates into T(hot blackbody) =T (cool blackbody), implying important changes of 
temperature for the same polar cap, indeed a small portion of the neutron star surface. 

4. Discussion 

The phase-resolved spectroscopy of PSR B0656+14, PSR B 1055-52 and Gcminga allows 
for a new view of the phenomenology of these middle-aged neutron stars. While the phase- 
averaged spectra of our targets look very similar, their phase-resolved behaviour is quite 
different. This is particularly evident from Fig. 12, Fig. 13, and Fig. 14, where the evolution 
of the cool and hot blackbody emitting radii, as well as of the power law flux, is shown as a 
function of the pulsars' phase. Owing to the not optimal quality of data collected for PSR 
B0656-I-14 and PSR B1055-52 above 2 keV, here we shall concentrate on the analysis of the 
two thermal components, presumably coming from the star surface, albeit from different 
regions at different temperatures. 

4.1. The cooler component 

First, let us examine the cooler component. It is seen to come from a region encom- 
passing a good fraction, if not the totality, of the neutron star surface as a result of the 
star cooling. We note that the best fitting emitting surface in the case of PSR B0656-I-14 
is very large (observed emitting radius of ~ 21 km) if compared with expectations for a 
standard neutron star. Distance uncertainties cannot ease the problem, owing to the very 
accurate radio VLBI parallax measurement (Brisken et al., 2003). In any case, a value of 
15 km, corresponding to a somewhat stiff equation of state (see Lattimer & Prakash, 2001 
and references therein), is marginally compatible with data (being allowed within the 99% 
confidence contour plot for and the emitting surface shown in Fig. 2). The pulsed frac- 
tion observed in the energy range where the cooler component dominates is of 14% for PSR 
B0656+14 (Fig. 6), of 16% for PSR B 1055-52 (Fig. 7), while for Geminga it is higher than 
30% (Fig. 8). 

We ascribe such a modulation to a variation of the emitting areas, going from a sizeable 
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fraction to the totality of the neutron star surface as a function of the rotational phase. A 
viable mechanism to provide some sort of phase- dependent "obscuration" of the bulk of the 
neutron star surface could be the magnetosphcric "blanket" , originally described by Halpern 
& Ruderman (1993) in their study of the soft thermal emission from Geminga (sec also 
Halpern & Wang, 1997): cyclotron resonance scattering by plasma in the magnetosphere at 
a few stellar radii could screen the thermally emitting surface during specific phase intervals, 
depending on the magnetic field configuration and viewing geometry. 

Alternatively, the flux variations could also be due to large-scale surface temperature 
modulations, expected as a consequence of anisotropic heat transfer from the neutron star 
interior (Greenstein & Hartke, 1983). As discussed in Sect. 3.3, such an interpretation is 
consistent with the observed spectral phase variation as well. With the current data we are 
not able to disentangle temperature and emitting surface variation without a complete phys- 
ical model. We note that the observed pulsed fraction values are much greater than expected 
on the basis of simple thermal models. As shown by Page et al. (1995), a few % modula- 
tion is expected in large-scale surface thermal emission, owing to the effects of gravitational 
bending. Pulsed fractions as high as the observed values could be explained assuming a 
peculiar beaming of the thermal surface emission. Owing to anisotropic radiative transfer 
in a magnetized plasma, an anisotropic angular emission pattern wrt the surface's normal is 
indeed expected, depending on temperature and magnetic field intensity and configuration^ 
(e.g. Zavlin & Pavlov 2002 and references therein). 



4.2. The hot spot(s) 

Turning now to the hotter component, we note that for all objects this is the most 
dramatically variable spectral component. It is natural to interpret such marked variations 
as an effect of the star rotation, which brings into view and hides one or more hot spots 
on the star surface. It is commonly accepted that neutron stars should have polar caps 
hotter than the rest of the surface. This could be due to different processes such as the 
bombardment of charged particles accelerated in the magnetosphere and faUing back to the 
polar caps along magnetic field lines (return currents, see e.g. Ruderman & Sutherland 1975; 
Arons & Scharlemann 1979), or anisotropic heat transfer from the neutron star core, which 
depends strongly on the magnetic field direction and is maximum along the magnetic field 
lines (see e.g. Greenstein and Hartke 1983). An association of the observed rotating hot 



^However, as noted in Sect. 3.1, current atmospheric models do not yield a satisfactory description of the 
spectra of the three musketeers 
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spot(s) with the neutron stars' polar caps seems therefore rather obvious. 

4-2.1. Luminosity and size 

The observed hot spot bolometric luminosities vary by more than 2 orders of magnitude, 
from the rather similar values found for PSR B0656+14 and PSR B1055-52 (5.7x10^^ erg 
s~^ and 1.6x10^^ erg s~^, respectively), to the much dimmer Geminga (1.6x10^^ erg s~^). 
The former two cases are in broad agreement with the theoretical expectations of Harding & 

Muslimov (2002) for polar cap heating due to downstreaming of e~^^~ generated by curvature 
radiation photons. In the case of the Geminga pulsar, close to the death line for creation 
of c"*"/" couples via curvature radiation, the lower polar cap luminosity is consistent with 
expectations for polar cap heating due to bombardment by particles created by inverse 
Compton scattered photons only (Harding & Muslimov, 2002; see Caraveo et al., 2004). 

Straight estimates of neutron star polar cap size, based on a simple "centered" dipole 
magnetic field geometry (polar cap radius Rpc = R\J where R is the neutron star 
radius, Q is the angular frequency and c is the speed of light), predict very similar radii for 
the three neutron stars, characterized by similar periods (233 m for PSR B0656+14, 326 m 
for PSR B 1055-52 and 297 m for Geminga, assuming a standard neutron star radius of 10 
km). The observed emitting radii are instead markedly different, with values ranging from 
~ 50 m for Geminga to ~2 km for PSR B0656+14. Psaltis et al.(2000) showed that the 
observed polar cap radius may be different from the actual one by a large amount, due to 
geometrical effects. Indeed, Caraveo et al.(2004) proposed a highly inclined viewing angle 
to reduce the surface of the emitting region in the case of Geminga. PSR B0656+14 and 
PSR B1055-52 face a completely different situation, since their polar caps are significantly 
larger than expected. Standard estimates, based on pure geometrical considerations, are 
clearly unsatisfactory. Even using the unrcalistically large neutron star radius found for PSR 
B0656+14, the inferred polar cap dimension is <600 m, far less than the hot spot radius. 
It is also possible, as suggested by Ruderman (2003), that thermal photons be significantly 
reprocessed higher up in the magnetosphere, interacting with charged particles. In such a 
picture, the same mechanism possibly providing the bulk of the cool blackbody modulation 
(via the phase-dependent screening quoted in Sect. 4.1) would obviously bias all emitting 
radii (and temperatures) measurements based on blackbody fits. 



4-2.2. Modulation with phase 



The observed phase modulation of the hot blackbody component is largely different for 
the three neutron stars both in amplitude and in pattern. In the case of PSR B1055-52 
a 100% modulation is seen, the hot spot disappears for 4/10 of the pulsar period. Such 
phenomenology strongly argues for an oblique rotator seen at high inclination. This would 
favour the historical Rankin (1993) interpretation of the radio polarization pattern of this 
pulsar (orthogonal rotator) versus the older one (almost ahgned rotator) of Lyne & Manch- 
ester (1988); however, the presence of only one visible "pole" would represent a challenge for 
a classical orthogonal rotator seen at high inclination. Conversely, the much lower modula- 
tion observed in the case of PSR B0656+14 seems consistent with an almost aligned rotator, 
with the polar cap always in sight. This picture is in good agreement with the radio pulse 
polarization results of Everett & Weisberg (2001), as well as with the above mentioned stud- 
ies of Rankin (1993) and Lyne & Manchester (1988). The remarkably high pulsed fraction 
of PSR B1055-52 (~70% in the 0.7-1.5 keV range) requires a significant beaming of the hot 
thermal component emission. Indeed, Psaltis et al.(2000) estimated that under standard 
assumptions about the star mass and equation of state, pulsed fraction higher than 35% 
cannot be produced even in the most "optimistic" case of an orthogonal rotator having very 
small polar caps with a high temperature contrast wrt. the rest of the star surface. The 
same consideration applies to the case of Geminga (pulsed fraction of ~ 55% in the 0.7-2 
keV range, but with a significant contribution of non-thermal emission) . In all cases a single 
peak per period is observed. This may suggest that we are seeing a single hot region on the 
surface of the star. The visibility of a single pole along the star rotation is consistent with the 
geometrical interpretation of the phase-resolved behaviour of PSR B0656-I-14, while for PSR 
B 1055-5 an important role of beaming should be invoked to explain why the opposite pole 
emission is not seen, in spite of the effects of gravitational light bending (Zavlin et al., 1995). 
An alternative possibility could be a magnetic field configuration different from a standard 
centered dipole. In the case of a magnetic field with significant multipole components, two 
(or more) hot polar caps could be very close on the star surface and their emission could be 
blended in a single peak. 



4.3. Phase alignement between the thermal components 

Owing to the sensitivity of heat transfer to the magnetic field, a temperature anisotropy 
should exist on the neutron star surface, with temperature increasing towards the magnetic 
poles (e.g. Greenstein & Hartke, 1983). Assuming the angular distribution of the radiation 
emitted by a surface element to be peaked along the normal to the surface and decreasing 
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towards larger angles (see Harding & Muslimov 1998 and references therein), we would expect 
the overall emission of the neutron star to be modulated as the poles come in and out of 
our line of sight. In particular, maximum flux should be observed when the line of sight is 
best aligned with the hotter regions. If the hot spots are indeed at the magnetic poles, a 
definite correlation between the hotter and colder components should be clearly visible in 
the phase evolution of the emitting regions. While this is indeed the case for PSR B1055- 
52, the contrary is true for PSR B0656+14, which shows a clear anti-correlation. Geminga 
shows an intermediate phenomenology, with a phase difference of ~0.25 between the hot 
and cool components (Fig. 14). Thus, one (or more) of the assumptions within the above 
simple scenario (temperature distribution resulting from a centered dipole magnetic field; 
pencil beaming) arc not correct. Different hypotheses should be considered: the surface 
temperature distribution is possibly more complicated, e.g. as a consequence of a multipolar 
magnetic field; the hot spots may not be located close to the center of the hotter surface 
region; the emission beaming function may have strong peaks at angles away from the normal 
to the surface element. We note that magnetospheric reprocessing of thermal photons (sec 
previous sections) could possibly ease the problem. Such a mechanism would introduce 
phase delays between thermal spectral components which would not be directly related to 
the properties of the surface, but driven by the magnetospheric structure. 

The observed phase-resolved behaviour of the three musketeers does not seem to argue 
in favour of the simple, "canonical" picture of neutron stars as inclined, rotating, centered 
dipoles. 

5. Conclusions 

In general, the new X-ray phenomenology revealed by the EPIC/XMM-Newton combi- 
nation presents new aspects of these three isolated, local neutron stars. Rotating, polar hot 
spots are clearly detected for the first time. It is tempting to link their origin to energetic 
particle bombardment, following the Ruderman mechanism, especially because certainly two 
(and possibly all) of the three musketeers are strong 7-ray sources. In the case of Geminga, 
moreover, Caraveo et al. (2003) have found a truly "smoking gun" evidence for this: the 
object's bow shock is filled with electrons of just the right energy (~ 10^^ eV) and luminosity 
(a few 10^^ erg s~^) to be the particles escaping from a polar cap. 

However, the measurement of polar cap sizes, rendered accurate by the excellent distance 
measurements, does not match simple theoretical expectations. 3-D geometry, beaming 
physics as well as the role of magnetospheric scattering will have to be invoked for detailed 
models. The same holds true for the apparent sizes of the cool blackbody components. At 
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least in one case they exceed the expected size of a standard neutron star. Simple atmospheric 
models worsen the problem. However, before considering the implications of this finding on 
neutron stars equations of state, we must understand if the parameters of the neutron stars 
surfaces can be reliably derived from their thermal radiation or if some other mechanism is 
biasing our results. 

Finally, the apparent puzzle posed by the difference in the cool and hot components 
phase alignements in the three objects is difficult to reconcile with existing pictures. This is 
rcminescent, however, of the situation at 7-ray energies, with different light curves. It will be 
very interesting to look for contemporary X-ray and 7-ray data when AGILE and GLAST 
will be in orbit with XMM-Newton. 
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Table 1: Journal of XMM-Newton observations. 



Pulsar /Date /Obs.Time 


Camera(mode)*^ 


Good Time 


Energy 


Photons(%bkg) 


Count rate 


PSR B0656+14 


pn(SW) 


5970 


0.15-8.0 


44600(1.7%) 


7.34±0.05 


2001-10-23 


pn(Ti) 


16850 


0.15-2.0 


120000(6.3%) 


6.67±0.04 


41.0 ksec 


MOSl(FF) 


37800 


0.15-8.0 


28100(2.1%) 


0.728±0.007 


PSR B1055-52 


pn(Ti) 


61900 


0.15-6.0 


84450(14.4%) 


1.167±0.009 


2000-12-14/15 


MOSl(FF) 


71000 


0.15-8.0 


17350(1.(3 ^/{) 


0.230±0.003 


81.4 ksec 


M0S2(FF) 


74250 


0.15-8.0 


18700(1.6%) 


0.247±0.003 


Gcminga 


pn(SW) 


55000 


0.15-8.0 


52850(5.4%) 


0.909±0.008 


2002-04-05 


MOSl(FF) 


76900 


0.15-8.0 


10170(2.2%) 


0.129±0.002 


103.3 ksec 


M0S2(FF) 


77400 


0.15-8.0 


11300(2.4%) 


0.142±0.002 



"^SW: Small Window; Ti: Fast Timing; FF: Full Frame 

Note. — Starting from the left, the column report (1) the target name, the date of the observation and 
the total time span (in ksec) of the observation; (2) the detector and its readout mode; (3) the good time 
(in sec) of the observation; (4) the energy range considered (in keV); (5) the overall number of counts in the 
source extraction region (see Sect. 2.2) and the fraction of background events in the specified energy range; 
(6) the background-subtracted count rate. 
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Table 2: Results of Timing analysis. 



Pulsar 


P Observed (ms) 


P expected (ms) 


PSR B0656+14 


384.9029(2) 


384.90300043(5) 


PSR B1055-52 


197.111812(5) 


197.111809432(8) 


Geminga 


237.1012(1) 


237.1012153(1) 



Note. — For each target, the best period, as computed from the EPIC X-ray data, is shown (column 2) 
together with the value expected on the basis of the extrapolation of published ephemeris: Kern et al.(2003) 
for PSR B0656+14, the ATNF data pulsar archive (http://www.atnf.csiro.au/research/pulsar/psr/archive/) 
for PSR B1055-52 and Jackson at al.(2002) for Geminga. The uncertainty quoted between parentheses refers 
to the least significant digit. 
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Table 3: Results of phase-integrated spectroscopy. 





PSR B0656+14 


PSR B1055-52 


Geminga*^ 


Njy (lO^o cm- 








4.3±0.2 


2.7±0.2 


1.07 (fixed) 


kTcBB (K) 








(6.5±0.1)xl05 


(7.9±0.3)xl05 


(5.0±0.1)xl05 


RcBB (km) 












8.6±1.0 


^Thbb (K) 








(1.25±0.03)xl06 


(1.79±0.06)xl06 


(1.9±0.3)xl06 


Rhbb (m) 








1800±150 


460±60 


40±10 


r 








2.1±0.3 


1.7±0.1 


1.7±0.1 


Ipz, (ph cm"^ 


s-^ keV 


-1 (0 


i 1 keV) 


4.3l?j X 10-5 


1.9+°:^ X 10-5 


6.7±0.7xl0-5 


Fo.2-8keV^ (erj 


5 cm~^ s" 






1.05x10-1^ 


2.2x10-12 


2.3x10-12 


Lpl'^ (erg s"^ 


) 






1.8x10^° 


8.1x10^0 


1.2x10^0 


Lcbb"^ (erg s~ 


') 






5.8x10^2 


4.4x10^2 


3.2x10^1 


Lhsb*" (erg s" 








5.7x10=^1 


1.6x10^1 


1.6x10^9 


Normp„^ 








Ti:l (fixed); SW:1.05 


1 (fixed) 




NormMOs/ 








0.96 


0.98 




NormM052^ 










1.07 




xVdof 








1.11 


1.02 


1.19 


dof 








368 


327 


73 



"Results from Caravco ct al.(2004) 
^Observed flux, 0.2-8 keV 

"^0.5- 10 kcV luminosity of power law eomponent 
''Bolometric luminosity of the Cool blackbody component 
^Bolometric luminosity of the Hot blackbody component 

■''Normalization factor to account for calibration differences; Ti: Fast Timing mode; SW: Small Window mode 

Note. — To compute luminosities we assumed a distance of 288 pc for PSR BG656+14 (Brisken et al., 
2003), of 750 pc for PSR B1055-52 (Kramer et al., 2003) and of 157 pc for Geminga (Caraveo et al., 1996) 
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Fig. 1. — Unfolded phase- integrated spectrum of Geminga. Only data from pn camera are 
plotted. This figure is adapted from Caraveo et al.(2004); a different color code is used here 
(see also Fig. 2 and Fig. 3). The best fitting spectral model is represented by the light blue 
line. As discussed in the text, this is obtained by the sum of a cool blackbody component 
(green), a hot blackbody component (red) and a power law (blue). Detailed values of the 
best fitting parameters are reported in Table 3. The inset shows confidence contours for the 
interstellar column density vs. the emitting surface for the cool blackbody. 68%, 90% 
and 99% confidence levels for two parameters of interest are plotted. Caraveo et al.(2004) 
fixed the N^^ value to 1.07x10^° cm~^ (resulting from ROSAT data) to obtain their set of 
best fitting parameters (see Table 3). 
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PSR B0656+14 




Energy (keV) 

Fig. 2. — Same as Fig. 1 for the case of PSR B0656+14. Data from pn (both Small Window 
and Fast Timing mode) and MOSl are plotted (black points). Detailed values of the best 
fitting parameters are reported in Table 3. 
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PSR B1 055-52 




1 2 
Energy (keV) 

Fig. 3. — Same as Fig. 1 for the case of PSR B1055-52. Data are from pn, MOSl and M0S2. 
See Table 3 for details on the best fitting spectral parameters. 
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Energy (keV) 



Fig. 4.— Unfolded spectra of PSR B0656+14 (red, data from pn and MOSl), PSR B1055-52 
(green, data from pn, MOSl and M0S2) and Geminga (black, data from pn). See text for 
details. 
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Fig. 5. — Light curves of the three musketeers. Data from pn (energy ranges: 0.15-2 keV for 
PSR B0656+14; 0.15-6 keV for PSR B1055-52; 0.15-8 keV for Geminga) have been folded 
using the radio timing solutions reported in Table 2 for PSR B0656-I-14 and PSR B1055-52 
and the EGRET 7-ray ephemeris (also reported in Table 2) for Geminga. The phase has 
been set in order to put the X-ray maximum at phase 0. The phases of the radio peaks have 
been marked with vertical dashed hues; their uncertainty is estimated to be ~0.01 (1/5 of 
phase bin) for the case of PSR B0656+14 and ~0.003 (~l/20 of phase bin) for PSR B1055- 
52. See text for further details. For PSR B1055-52 "radio peak 1" refers to the highest 
peak in the radio profile, see also Fig. 7. We plotted also the phases of the 7-ray peaks for 
Geminga; however, as discussed by Caraveo et al.(2004), the propagation of errors in the 
extrapolation of the EGRET ephemeris makes their position uncertain by ±0.15 ( such a 
phase interval corresponds to the gray-shaded regions; "7 peak 1" refers to the highest peak 
in the 7-ray profile, see also Fig. 8) 
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PSR B065B+1+ 




0.15-0.50 keV 
Cool blaekbody 
PF=(14.3±0.5)% 

0.50-0.70 keV 
Intermediate range 
PF=(4.5±1.0)% 

0.70-1.50 keV 
Hot blockbody 
PF=(21.3±2.5)% 

1.50-7.00 keV 

Power law 
PF=(75±20)% 

Radio 0.S7 GHz 



0.5 1 1.5 

Relative phase 



Fig. 6. — Lightcurves of PSR B0656+14 in different energy ranges. Data (obtained with 
pn Fast Timing mode observations, with the exception of the hardest band, based on pn 
Small Window mode data) have been folded using the radio ephemeris quoted in Table 2. 
The alignement in phase is the same chosen for Fig. 5. Although always single peaked, the 
pulse profile changes significantly going from the softest energy range (dominated by the cool 
blackbody) to the 0.7-1.5 keV range (dominated by the hot blackbody). Note the minimum 
in the pulsed fraction in the intermediate 0.5-0.7 keV range. Above 1.5 keV the lack of 
statistic hampers a detailed study of the pulse profile. The radio hght curve is shown in the 
lower panel. The uncertainty on the phase alignment of the X-ray light curve with the radio 
one is of ~0.01. The single radio pulse is seen to trail the X-ray maximum by ~ 0.2 in phase. 
The phases of the two optical peaks (Kern et al., 2003) as well as of the 7-ray peak claimed 
by Ramanamurthy et al.(1996) are also shown. See text for further details. 
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PSR Bl 055-52 



0.15-0.35 keV 
Cool blackbodji 
PF=(16.0±1.0)% 



0.35-0.70 keV 
Cool blackbody 
PF=(29.0±1.0)% 

0.70-1.50 keV 
Hot blackbody 
PF=C67±3)% 



1.50-6,00 keV 
Power law 

pr=(90±io)% 



Radio 0.67 GHz 



0.5 1 1.5 2 

Relative phase 

Fig. 7. — Light curves of PSR B1055-52. Note the broader profile of the X-ray pulse below 
0.35 keV and the small (~ 0.1) phase shift wrt. higher energies. Note the value of the 
pulsed fraction, which grows with energy and is remarkably high (~ 70%) in the 0.7-1.5 

keV range, where the neutron star emission is dominated by the hot blackbody component. 
For a more detailed study of the energy resolved pulsed fractions see Becker & Aschenbach 
(2002). Above 1.5 kcV the low signal-to-noise hampers a detailed study of the pulse profile. 
Wc plotted in the lower panel the radio light curve. The uncertainty on the phase alignment 
is of ~ 0.003. The two radio peaks are observed to bracket the single X-ray peak. The 
7-ray pulse (Thompson et al., 1999) occurs in the gray-shaded phase interval between the 
two vertical dashed lines. 
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Fig. 8. — Light curves of Geminga. As discussed by Caraveo et al.(2004), the pulse shape 
chang function of energy. The single, broad peak observed at low energy (where 

emission from the cool blackbody dominates) changes to two peaks at higher energies (where 
the power law component dominates). The pulsed fraction is maximum in the 0.7-2 keV 
range, where the hot blackbody component is more important. The lower panel show the 
EGRET 7-ray lightcurve. The extrapolation of the 7-ray ephemeris makes the absolute 
phase ahgnement uncertain by ±0.15, owing to the long time span between the EGRET and 
the EPIC observations. Note that the numbering of the phase intervals defined in the upper 
panel is different from that used by Caraveo et al. (2004). Their phase 1 is here phase 3 and 
similarly for all. 
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Fig. 9a. — Phase-resolved spectra of PSR B0656+14. Photons have been selected in corre- 
spondence of the phase intervals marked in Fig. 6. Phase intervals from 1 to 5 are shown 
here; phase intervals 6-10 are shown in Fig. 9b. Upper panels on the left column present, 
for each phase interval, the observed spectrum (data points) compared to the best fit model 
(solid line) of the phase-integrated spectrum (upper plots); lower panels show the difference 
between data and such model in units of statistical errors. To ease the visibility of the de- 
viations of phase-resolved spectra from the averaged spectrum template, spectral channels 
have been rebinned. Panels on the right column show, for each phase interval, the unfolded 
spectrum together with its best fit model. The model components are also plotted (color 
code as in Fig. 2). An animated version of Figures 9a/b, lOa/b and lla/b is available at 
http: / / www.mi.iasf.cnr.it/~deluca/3musk/ 
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Fig. 9b. — Phase-resolved spectra of PSR B0656+14. Same as Fig. 9a, phase intervals 6-10 
are shown. 




0.5 1 2 5 0,5 1 2 5 

Energy (keV) Entrgy (keV} 



Fig. 10a. — Phase-resolved spectra of PSR B 1055-52. Phase intervals 1-5 (according to the 
notation of Fig. 7) are displayed. See caption to Fig. 9a for explanations. Data have not been 
rebinned. Note the much higher deviations of the phase- resolved spectra wrt. the average 
one. 
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Fig. 10b.— 



Same as Fig. 10a, phase intervals 6-10 are shown. 
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Fig. 11a. — Phase- resolved spectra of Geminga. Phase intervals 1-5 (see Fig. 8) are shown 
here. See caption to Fig. 9a for explanations. The figure is adapted from Caraveo et al.(2004), 
according to the phase numbering and color code adopted in this work. 
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Fig. lib. — Same as Fig. 11a, phase intervals 6-10 are displayed. 
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Fig. 12. — The parameters best fitting the phase-resolved spectra of PSR B0656+14 are 
plotted as a function of the pulsar phase defined as in Fig. 5. Both cool and hot blackbody 
emitting surfaces evolve throughout the pulsar phase following a sinusoidal pattern showing 
an overall ~10% modulation (wrt. the average values) on the emitting radii value. Note 
the marked anti-correlation between panel 1 and panel 2 with the hot blackbody peaking in 
correspondence of the cool blackbody minimum. The power law component has a different 
phase trend wrt. the thermal components, with a single, sharper peak. 
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Fig. 13. — Same as Fig. 12, for PSR B1055-52. While the cool blackbody emitting radius 
shows a <10% modulation (wrt. the average value), the hot blackbody component show 
a dramatic, 100% modulation since its contribution is null in 4 of the 10 phase intervals. 
Note that for PSR B1055-52 the two thermal components have a similar time evolution, 
with a phase shift as low as ~0.1. The power law component has a different modulation; it 
is difficult to assess wether its profile is single-peaked or double-peaked, owing to the lower 
signal to noise in the high energy portion of the spectra. 
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Fig. 14. — Same as Fig. 12 and Fig. 13 for the case of Geminga. The figure is adapted from 
Caraveo et ah (2004) (see their Figure 4), according to the choice of phase adopted in this 
work. The cool blackbody component shows a ~ 15% modulation (wrt. the average value) 
of its emitting radius, with a sinusoidal profile. Conversely, the hot blackbody is 100% 
modulated, and disappears for 1/10 of the pulsar period. The power law component has a 
remarkably different, double-peaked phase profile and shows a significant unpulsed fraction. 



